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^ . Abstract 
Ch 

The aim of this paper is to tackle the nonlinear optical reconstruction problem. 
Given a set of acousto-optic measurements, we develop a mathematical framework 
£ — ■ for the reconstruction problem in the case where the optical absorption distribution 

is supposed to be a perturbation of a piecewise constant function. Analyzing the 
acousto-optic measurements, we establish a new equation in the sense of distributions 
for the optical absorption coefficient. For doing so, we introduce a weak Helmholtz 
decomposition and interpret in a weak sense the cross-correlation measurements using 
the spherical Radon transform. We next show how to find an initial guess for the 
unknown coefficient and finally construct the true coefficient by providing a Landweber 
type iteration and proving that the resulting sequence converges to the solution of the 
system constituted by the optical diffusion equation and the new equation mentioned 
above. Our results in this paper generalize the acousto-optic process proposed in [¥] 
for piecewise smooth optical absorption distributions. 
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1 Introduction 

Let be a bounded C 1 -domain of ]R d , where d E {2,3}. We denote by v the outward 
S_i . normal to d£l, the boundary of f2. We need the following functional spaces. For m a 

non- negative integer, we define the space H m (Q) as the family of all functions in L 2 (fi), 
whose weak derivatives of orders up to m also belong to L 2 {VL). For m > 1, the space 
£P»-i/2(an) denotes the set of the traces on <9s2 of all functions in H m (Q). We let flj*(fi) 
be the closure of C£°(0) in H m (£l), where C£°(fi) is the set of all infinitely differentiable 
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functions with compact supports in f2. We denote by H m (Q) the dual of H™(Q). Finally, 
for p > 1, we introduce W m,p {0) as the space of functions whose weak derivatives of orders 
up to m are functions in L p (n) and W™' p (n) to be the closure of C£°(ft) in W m ' p {n). Here, 
L p (n) is defined in the usual way Note that W m > 2 {ti) = H m (Q) and W™' 2 (Q) = H^(Q). 

Suppose that Q represents an optical medium and let a* : Q — > M + be the optical 
absorption coefficient of f2. When the medium Q is illuminated with infrared light spots, 
the optical energy density € H 2 (Q) inside $1 satisfies the diffusion equation 

-A$* + a*<£>* = in n, 

(1.1) 

ld v <&* + ^ if = g on dQ, 

where I > is the extrapolation length, computed from the radiative transport theory 
[25], and the illumination function on the boundary g 6 H l l 2 (dSl) satisfies g > a.e. on 
d£l, and d u denotes the normal derivative at dVl. 

In diffuse optical tomography, the inverse problem is to reconstruct the optical absorp- 
tion distribution a* from measurements of the outgoing light intensity on d£l given by 
9u^*\dCit see [Hi [27]. It is worth mentioning that, in our diffusion equation model (jl.ip . if 
I 0, then knowing <1>* or d v Q* on dVt is mathematically the same. 

Diffuse optical tomography produces images with poor accuracy and spatial resolution. 
It is known to be ill-posed due to the fact that the outgoing light intensities are not very 
sensitive to local changes of the optical absorption distribution j9[ [191 EH EZ] ■ In [5] we 
have proposed an original method for reconstructing the optical absorption coefficient by 
using mechanical perturbations of the medium. While taking optical measurements the 
medium is perturbed by a propagating acoustic wave. Then cross-correlations between 
the boundary values of the optical energy density in the medium changed by the prop- 
agation of the acoustic wave and those of the optical energy density in the unperturbed 
one are computed. Finally, under the Born approximation [12] . the use of a spherical 
Radon transform inversion yields a reconstructed image for a*, which has a resolution of 
order the width of the wave front of the acoustic wave propagating in the medium. The 
Born approximation linearizes the reconstruction problem. It consists of assuming that 
a* is close to a constant and taking the background solution of the diffusion equation for 
constant optical absorption in place of <£* [27J as the driving optical energy density at 
each point in O. 

The idea of mechanically perturbing the medium has been first introduced in [1] for 
electromagnetic imaging. On the other hand, it is also worth emphasizing that this ap- 
proach is different from the imaging by controlled perturbations [21 [31 [131 El El ESI [29] . 
where local changes of the parameters of the medium are produced by focalizing an ultra- 
sound beam. Both techniques lead to resolution enhancements. In imaging by controlled 
perturbations, the resolution is of order the size of the focal spot while here it is of the 
order of the width of the wave front of the wave propagating in the medium. 

This paper aims to generalize the acousto-optic process behind the Born approxima- 
tion. We tackle the nonlinear optical reconstruction problem. We develop a mathematical 
framework for the reconstruction problem in the case where the optical absorption dis- 
tribution is a perturbation of a piecewise constant function. We introduce an iterative 
reconstructing algorithm of Landweber-type and prove its convergence and stability. For 
doing so, we introduce a weak Helmholtz decomposition and interpret in a weak sense the 
cross-correlation measurements. 
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To describe our approach, we employ several notations. Each smooth component of o* 
is called an inclusion. The background of a* is assumed to be a known positive constant 
and denoted by ao- Assume further the knowledge of a lower bound a and an upper bound 
a of a t , both of which are positive. Finally, let D <s Q be known and such that 

a* = ao in Q \ D. (1-2) 

We next impose some conditions on the unknown inclusions. Let k > 1 denote the number 
of inclusions and Ai be occupied by the ith inclusion. Assume: 

I\. for any i £ {1, . . . , k}, Ai is a smooth subdomain of 0, dAi is connected; 

li- for any j ^ i, Aid Aj = 0; 

h. \J k i=l Ai cs D. 

All of the assumptions above suggest the definition of the class (A) C L°°(il), which 
contains a*. 

Definition 1.1 The function a is said to belong to class (A) iff there exist k > 1, 
Ai, ■ ■ ■ , (s D satisfying I\,l2 and I3 and a\, ■ ■ ■ , € C 2 (Ai, [a, a]) such that 

k 

a = '^2a i l Al , (1-3) 

where, again, ao was introduced in Aq = $7 \ uf =1 Aj ; and 1a x denotes the charac- 

teristic function of Ai . 

Our main results in this paper can be summarized as follows. A spherical acoustic 
wave is generated at y outside Q. Its propagation inside the medium Q changes the 
optical absorption distribution. Due to the acoustic wave, any point x € moves to its 
new position x + Vy, r (x), where Vy >r is defined by f)3. 1 j) with r being the radius of the 
spherical wave impulsion. By linearization, the displacement field is approximately Vy )T as 
the thickness rj of the acoustic wavefront goes to zero. Hence, the optical absorption of the 
medium changed by the propagation of the acoustic wave is approximately a*(x + Vy, r ), 
up to an error of order 77. 

Using cross-correlations between the outgoing light intensities in the medium changed 
by the propagation of the acoustic wave and those of in the unperturbed one, we get the 
data M v (y, r) given by f)3.5[) . In Propositions 14.21 and !4. 31 we show that M rj (y, r) converges 
in the sense of distributions to M(y, r) as rj — > 0. We refer to M(y, r) as the ideal data. 
Making use of a weak Helmholtz decomposition, stated in Lemma [2,61 we relate in Theorem 
15.11 the ideal data to the gradient of ^Va*. Since a* is piecewise smooth, Va* can be 
defined only in the sense of distributions. Technical arguments and quite delicate estimates 
are needed in order to establish the fact that the gradient part of ^Va* can be obtained 
from the cross-correlation measurements using the inverse spherical Radon transform. 
Based on this, we propose an optimal control approach for reconstructing the values of a* 
inside the inclusions. For doing so, we first detect the support of a* — 00 as the support of 
the gradient part of the data ^Va*. In fact, Lemma [231 shows that the support of the 
data yields the support of the inclusions. Their boundaries are detected as the support 
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of the discontinuities in the data. Proposition 16,11 provides a Lipschitz stability result for 
reconstructing piecewise constant optical absorption. In contrast with the recent results in 
[U llCH fTT]. Proposition ^ . 1 1 uses only one measurement but the supports of the inclusions are 
known. Minimizing the discrepancy functional (j6.3|) we obtain the background constant 
values of the optical absorption inside the inclusions. Next, in order to recover spatial 
variations of a* inside the inclusions, we minimize the discrepancy between the linear 
forms F[a] and Ai/j given by (|6.6j) and (|6.13p . respectively. We prove in Theorem 16.31 
that the Frechet derivative of the nonlinear discrepancy functional is well-defined and 
establish useful estimates as well. We introduce an iterative scheme of Landweber-type for 
minimizing the discrepancy functional and prove in Theorem 16.51 its convergence provided 
that the optical absorption coefficient is in the set K defined by (|6.4p . 

2 Preliminaries 

2.1 Some basic properties 

We first recall the following results. 

Proposition 2.1 (weak comparison principle |5j) Let a € L°°(Q) be a nonnegative 
function and assume that $ £ H (fl) satisfies 

( + a $ > inn, 

| + > ondn. 

We have <& > a.e. in VI. 

Lemma 2.2 (Lemma 4.1 in |5j) Let D be as in and assume that g € H l l 2 {d£l) 

is nonnegative. There exist two positive constants A and A such that for all a € (-4), the 
solution $ of 

-A$ + a$ = in O, 

(2.2) 

ld u $ + $ = g on dtt, 
satisfies 

A < $ < A in D. (2.3) 



Lemma 2.3 (Lemma 4.2 in |5j) Let T be the map that sends a € (-4) into the unique 
solution of with a replacing a*. Then, T is Frechet differentiable. Its derivative at 
a is given by 

DT[a](h) = <p, (2.4) 

for h € L°°(p,), where (p solves 

I — A09 + aip = —hT\a] in Q,, 

^ (2.5) 
ld u ip + ip = on dQ. 
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Moreover, DT[a] can be continuously extended to L 2 (Q) by the same formula given in 
(gg j and (KW w tth 

\\DT[a}\\ c{L 2 {nlH i m <CA, (2.6) 

where A is defined in Lemma \2.2\ and £(L 2 (Q) , H 1 (£1)) is the set of bounded linear opera- 
tors from L 2 (Q) into H 1 (Q). 

The following lemma will be helpful to prove the uniqueness of the constructed coeffi- 
cient. We refer to Appendix [X] for its proof. 

Lemma 2.4 Let £1' be the union of several subdomains of Q such that fi \ f2' is path 
connected. If (j) is a bounded solution to 

[ -A6 + cd> = inQ\n', 

<^ (2.7) 
I ld u <fi + 4> = on dfl, 

for some nonnegative constant c and d u (f> = on d£l, then = in Q \ f2'. 

Corollary 2.5 Let Aq, A±, ■■■ , Ak be as in Definition li.il and let a G (A) be defined by 
such sets. Denote by cpj, j = 1, • ■ • , k, the solution of 

{— Aifj + a(fj = lAj& in 
ld u (fj + if j = on <9f2, 

with $ being the solution of H2.2\) . Then, the set {d u (pj\gn} is linearly independent. 
Proof. Define 

k 

if = a jfj m 

for some «i, • • • , a& E R, and assume that d u ip = on d£l. It is obvious that (p is the 
solution of 

k 

—Aip-\-a<p = a j ^-Aj ^ in O, 

3=1 

ld v >*p + (p = on <9f2, 

and, hence, satisfies (|2.7p with c = ao and 0' = U^ =1 Ai. Thus, by Lemma l2~il tp = in 
fit). On the other hand, for each i € {1, • • ■ ,k}, <p solves 

j -Aip + aiip = ailAi® in A i} 
} <p = on dAi. 

We can now apply the strong comparison principle (see, for instance, Lemma 3.1 in |22j) 
and the Hopf lemma to see that d u (p ^ on dAi. This contradicts to the fact that = 
in Aq. □ 
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2.2 The Helmholtz decomposition in the sense of distributions 

The Helmholtz decomposition plays a crucial role in when we established a differen- 
tial coupling system for a, where a was supposed to be in C 2 (0). Fortunately, when a is 
no longer smooth but $ 2 Va belongs to C H~ 1 (£l) d for all € C 1 ($7), a corre- 

sponding Helmholtz decomposition remains true. Note that for all a G (.4.) and <E> G C 1 (J7), 
$ 2 Va G (H 1 ^)* in the sense that 

($ 2 Va,f) (H i( n) d ) , )jH -i (n )d = ($ 2 V(a - a ),v)(Hi («)<»)* ,/ri(n)<* 

= - / (a - a )V • ($ 2 w) dx. (2.8) 
./£> 

The domain of the integral above is written as D instead of O because a — ao = in O \ Z), 
where D is introduced in (|1.2|) . By the same reason, we do not require the boundary zero 
value for the admissible test functions. The last equation in (|2.8p suggests that it might 
be sufficient to impose a G C 1 (^4j), instead of C 2 (Ai), i = 1, • • • , k, as in Definition 11.11 
However, we need the differentiability of a up to second order in each inclusion for some 
later regularity and estimation purposes. 
The following result holds. 

Lemma 2.6 For any U in H~ l (Q) d there exist tp G L 2 (0) and \& G H~ l (ft) d such that 

U = W + * 

with V • $ = 0. in particular, if U = <£ 2 Va /or some a £ i i/ien "0 *s continuous and 
discontinuous at the point where a is, respectively. 



Proof. Letting U = (U\, • • • , Ud) G H l {Vl) d , we denote by u = («i, • • • , ii<i) the solution 
of 

-An = U in Q, (2 9) 



it 



on 90. 



The vector u € i?o(^) d is actually the Riesz representation of U in i7g($7) d . Applying the 
classical Helmholtz decomposition for u (see, for instance, |14j). we can find / G -ff^O) 
and G G ff(curl, O) := {w G L 2 (0) d : V x w G L 2 (0) d } such that 

« = V/|VxG. (2.10) 

Here, V • G = inside O and 

G x v = on 90. (2.11) 

Moreover, / is a solution of 

f A/ = V • u in O, 

< 2.12 
\ cV = on 90. 

Since u belongs to ff 1 (0) d , V • u G L 2 (Q). By standard regularity results, we see that 
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In view of (j2.9|) . taking the Laplacian of (|2.10p yields 



U = V^ + ^, 

in the sense of distributions, where ip = Af G L 2 (Q) and ^ is divergence free. 

We next prove the second statement of the lemma in which U = <3? 2 Va for some a G 
(A). The main tools we use here are the H 2 - and (^-regularity results. Fix j G {1, • • • ,d} 
and i G {0, • • • , k}. Denote by Uj the jth component of the vector u, defined in ()2.9|) . 
Since Uj G Hq(Q), it belongs to H (A{). The function Uj solves 

- A Uj = <f> 2 d Xj a, (2.13) 

in Aj. Applying Theorem 8.8 in [1 tj| . we see that tij is in .£f 2 (^) for all A\ <s -Aj. Hence, 
differentiating (|2. 13|) gives 

-Ad XlUj =d Xl ($ 2 d Xj a) 

in A' { for alH = 1, • • • ,d. Since d Xl Uj G -ff 1 ^-) and d Xl ($ 2 d Xj a) G L 2 ^), we can apply 
the (^-regularity result in [2D] to see that is in C l (A") for all A" d ^4.. This implies 

Uj G C 2 (Ai). Considering the differential equation in (|2.12p in each inclusion and following 
the same regularity process, we see that / G C 2 (Ai). Hence tp = Af is continuous in Ai, 
which is also the set of continuous points of a. On the other hand, since U = < I )2 Va 
involves Dirac distributions supported in UjcMj, V • u is not continuous across UjcMj, so 
are / and ^ = Af. □ 



3 The set of data 

In this section, we describe the set of data obtained by the acousto-optic process introduced 
in [3]. The basic idea in order to achieve a resolution enhancement in imaging the optical 
absorption distribution is a s follows. We generate a spherical acoustic wave inside the 
medium. The propagation of the acoustic wave changes the absorption parameter of the 
medium. During the propagation of the wave we measure the light intensity on d£l. The 
aim is now to reconstruct the optical absorption coefficient from such set of measurements. 

Let a G (A) represent the true coefficient a*. Let be the unit sphere in R rf . Let 
fi > and let = [iS x , the sphere of radius /i and center 0, be such that Q stays 
inside Su. We perturb the optical domain by spherical acoustic waves generated at 
point sources y G Sn. Let r G [ro,i?] be the radius of the spherical wave impulsion, where 
ro and R are the minimum and maximum radii so that the spherical waves generated at 
point sources on can intersect 0. Let r] <C 1 be the acoustic impulsion typical length 
representing the thickness of the wavefront. Let the position function P be defined by 

P '. x i y x ~\~ v^ r (x), x G f2, 

where 




and w is a smooth function supported on [—1,1] with ||itf||oo = 1- Here, || ||oo denotes 

II IU~a-i,i[)- 
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In [3], we have shown that the displacement function at the point x caused by the 
short diverging spherical acoustic wave generated at y is given by 

ul r {x) = p- 1 (x)-x, (3.2) 

Let C be the cylinder x [ro, R]. For each (y, r) € C, a u v r (x) denotes a(x + Uy >r (x)) 
and $ u j7 ^ is the optical energy density in the displaced medium, which satisfies 

f - A$,.n + o,.» <J> i =0 in fi, 

Physically, the outgoing light intensities B u ^\qq and d u <& u v Jqq are measured. We are thus 
able to assume the knowledge of the cross-correlation measurements: 

4 / g(du$ - d v $ u n r )da, yeS^,r>0. (3.4) 
T Jon 

Integration by parts shows that the quantity above is equal to 

M v(y> r ) = ^2 J ( a < r - a)$$ u v r dx, (3.5) 

which is considered as our set of data. Here, the coefficient l/rf is put in front of the 
integral because both and $ u ^ r are bounded (Lemma I2.2[) and 

ll«< r - a||£i ( fi) = 0(rj 2 ) as??^0 + , (3.6) 

provided that the following technical condition, named as (H), is imposed: there exists 
5 > such that for all x € dA{ fl T, v (y, r), either 

Hi : the angle formed by the ray x — y and the normal outward vector of Ai at x is greater 
than 5; or, 

H2: the curvature of dAi is different to that of the circle or sphere {z € M. d : \z — y\ = 
\x — y\} at x if the angle above is smaller than 8. 

Here, 

S^(y, r) = {z € M d : r — rj < \z — y\ < r + rj}. 
In fact, this condition guarantees that 

lAAP-^Ai)] + \AiAP(Ai)\ < 0( V 2 ). (3.7) 

Denote 

V £ (S) = {xeR d , 3yeS, \x-y\ < e}, (3.8) 

for any smooth surface S of M. d , and e > 0. Since 5 is smooth, the volume of V e (S) is given 
by 

V £ {S) = 2a{S)e + 0{e 2 ). 
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Fix (y, rf) G C and write 

n 

hu^r -0||l1(0) = J2 \ a u^r~ a \ dX 

= J2 / \ a u?, r - a \dx+ \a u v - a\dx. (3.9) 

As is supported on S^(y,r) and ||uy jr ||oo — ?7> 

/ l«^r~ a l da; = / |a< r -«l^ 

JAinP- 1 ^) Js^nAinP- 1 ^) 



< 



|Va i || L oo (Ai) o-(S'(0,.R))77 2 , 



where cr(5(0, i?)) is the surface measure of the sphere of center O and radius R. The 
second integral in (|3.9|) is bounded by 0(n 2 ) because of fj3.7|) and the boundedness of a. 



4 The behavior of M v as r] approaches + and the ideal mea- 
surements 

Consider the open cylinder C := S^x (0,R) with its classical product topology. 
The construction of 

M V :C K 
0/> r ) ^ ^2 J^ a <r- a )^ul r dx, 

has been described in this previous section. The knowledge of this function is obtained 
from those of <?, d v Q and d v § u v r on dVt. In this section, we study the limit of M v as 
r] — >■ + . This, together with a weak version of Helmholtz decomposition and the spherical 
Radon transform, will help us to detect all inclusions. 

Lemma 4.1 For any n > 0, M r] is a continuous map on C. 



Proof. It is sufficient to consider only the case r > ro because M v (y, r) = for all r < ro 
and y G S^. Fix (y,r) G x (r ,i?) and let {(y n , r„)} n >i C 5 M x (r ,-R) converge to 
(y, r). Noting that a« y r is continuous except on the zero measured set 

{x + u^ r (x) : x G U" =1 dA;}, 

we have 

a ( X + U y n ,r n ( x )) -> a O + U&rO")) 

a.e. in f2. On the other hand, since a is uniformly bounded, so is 

\ a ( x + U y n ,r n i X )) - a ( X + u yA X ))\ 2 - 
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It follows by the Lebesgue dominated convergence theorem that 

a..n — > a..ri in L 2 (Q.) 

as n — > oo. This implies 

<f>,.v -4 

a Vn,r n a y,r 

in both H 1 ^) and L 4 (f2). Note that the L 4 convergence above is valid because d is either 
2 or 3. A direct calculation yields 

\n 2 (M v (y n ,r n ) -M v (y,r))\ 

/ [(a,.* - a)$$..v - (a,.v - a)$$..v ]dx 

Jq a Vn,r n ' a Vn,rn v a y,r ' "y,r J 

-a|$|<I><? - <t>,.i \dx+ I \a,.v - a u n \®$,.v dx, 



< / a 



< 2a||$|| L 4 (n) ||$ u 2 nij . n -$ u3ir || L 4 (n) 



+ <V 



°n^r Hi 2 (^) \\®\\L±(Q) II ll^.r H J< 4 («) 



The lemma follows. 



□ 



Lemma 14.11 guarantees that is measurable. In the case that a is smooth, which 
has been studied in |H |5], M v (y, r) « Va ■ u^ r $ 2 when r\ is small. However, when a is 
piecewise smooth, we need to establish a similar approximation in the weak sense. The 
following proposition holds. We refer to Appendix [B] for its proof. 

Proposition 4.2 Let C = S^x (0, R). For any 0<rj<l, define the continuous function 

1 



M v (y,r)- 

T Jn 



(a - a )V • ($ 2 i# r )ds, (y, r) G C, (4.1) 
where $ is i/ie solution of with a replacing a*. Assume (H) holds and, consequently, 



3.1) is valid. Then there exists c > 0, independent of (y,r), such that 

M v (y,r)-M n (y,r) < cry, V(j/,r)eC7. 
It follows from Proposition 14.21 that for each (y, r) € C, 

lim M v (y,r)= lim M v (y,r) := M y>r . 

r/^0+ r;->0+ 



(4.2) 



(4.3) 



We cannot expect that M is a smooth function on C because Uy jr /rp, and hence Vy^/rj 2 , 
converges to a distribution supported on the circle (or sphere) S(y,r) = {z : \z — y\ = r}. 
The limit in (j4.3|) is understood as follows. 
Let 

G(C) = {f£L 2 (C):d r feL 2 (C)}, 
be a Hilbert space, endowed with the norm 



\G(C) 



+ 1 1 a. 



Il2(c)- 
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Let 7 be the (continuous) trace operator from C to S M x {0, R} and denote 

G (G) = 7 - 1 (0) = {/ G G(G) : 7 (/) = 0}, G^G) = G (G)*. 

We have the following relations 

H]{C) C G (C) C L 2 (G), L 2 (C) C G^G) C ir^G). 

Let || ||i denote || ||_Li(i— i in- The following is the main result of this section. It is a 
direct consequence of Proposition 14.21 

Proposition 4.3 The function converges to the ideal measurements M in G _1 (G) as 
rj — > + with 



{M, v)(c^(c)y,c^{C) 



= -ro\\w\\i 

where 



R 



d 

a{y + r£)- 



(9 



Proof. For any in Go(G), we have 



(m v ,<p) = -[ ! I \a-a Q ){x)V- x U\x) V Mt L v{y,Adxdrdy 

* t#, r (s) 



(a - a )(s) V $ 2 (x) / y ' 3 r)dr ) cfedy. 

2/€S M Jf2 V Jr=0 V 



Then by the change of variables x = y + p£ we get 

«y,rO&) _ r ,p-r 



w{- )t- 

rj A rrf n 



Hence we can write 

M v ,ip i 



[ [ l R (a-a )ty + P 0-?-( P d - 1 $ 2 (y + pO f ( ^ <p(y,r)dr) dpd^dy 

Js» Js d -i J P =o dp V Jr=o rv V V J ) 



'5 M is 1 *- 1 ip=0 
Since 



we deduce that 



1 p-r\ n^o .. .. 
-u; — > \\w \\i0p, 

V V V J 



r o (P~ r \ ( u v^o \\w\h r o , x 

— w[ \ip(y,r)dr — > <p{y,P), 

=0 rij V V J P 



and then 



»7-i>0 



w\Ur 



[ [ [ a(y + pO^-(p d - 2 ^ 2 {y + POv(y,P))dpdCdy, 
Js u Jn=0 Js"- 1 d P v J 



as desired. □ 
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5 Detecting the inclusions 

Using the fact that <£ 2 Va G C H~ 1 (Q) d , we can employ Lemma ETUI to write 

that 

$ 2 Va = W + *, (5.1) 

where \& is a divergence free field and ip € L 2 (il). Since that both <£ 2 Va and are in 
(H 1 (£l) d )* , so is ^. Moreover, it follows from the usual integration by parts formula and 
the boundary condition ()2. llj) that 

(tf,Vw) = 0, VueC°°(n). (5.2) 

For a distribution / € (Qj°(C))*, we define its spherical Radon transform TZ[f] in the 
sense of distributions by 

(^[/]>¥')(Cg o (C))* ) C§ o (C0 = (/,^*[^])(C^(C))*,eg°(C0, 

where 

K*[<p](x) = [ <p(y,\x-y\)dy for ^ G C °°(C7). 
We have the following result. 

Theorem 5.1 The spherical Radon transform 7Z[ip] of tp satisfies the equation 

M = r \H\ ir ^^ (5.3) 
or 



in the sense of distributions. 

Proof. Let ip G Cq°(C), and for a fixed y £ we define 

F 1/ (a;) = (p(y, \x-y\)- ^ x G O. 

| a; — y| z 

For any y G 5^, the vector F y is in H 1 ^) 11 because \x — y\ > r^. Equation (|5.ip yields 

(<S> 2 Va,F y ) = {ViP,F y ) + {*,F y ) 
and since i 7 ^ is the gradient of the function given by 

dp, 

P 

it follows from (|5.2p that (^,F y ) = 0. Here, ( , ) denotes the duality pair between H 1 (ft) d 
and (il 1 ^)*. Therefore, 

($ 2 Va,F 2/ ) = {V^,F y ). (5.4) 



A simple calculation shows 

($ 2 Va,F y ) = f aV-($ 2 F y )dx 
Jn 



R r 

/ [aV • ($ 2 Fj,)] (y + r^r^^dr, 
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and hence, 



<$ 2 Va,F y > 



R 



o Js d - 1 nfiy, r 



d r i 
<*(y + rO*- <$> 2 (y + r£)<{>(y,r)r d - 2 d^dr. (5.5) 



dr 



Combining (|4.4p . (|5.4p . and (|5.5[) implies 



(M,^) = ||w|| L i ( n) / (Vip,Fy)dy 



\ w \\LHn) / / V'V • (Fy)dxdy 
ls„ Jn 

R 



\w\\ L i (n) I I I [iPV ■ (F y )} (y + r£)r a 'd^drdy 

iSp Jo Js d - 1 nn Vtr 



ls» Jo Js^nn y , r dr L 

-fo || w || .z,i(n) / / K[^}(y,r)— Lp{y,r)r d ~' A drdy 

' Su, J 



d-2 



d^drdy 



= r \\w\\ L i {n) (r a 
and the proof is complete. 



dr 



,<P), 



□ 



Remark 5.2 Theorem \5.1\ provides the knowledge of the derivative of the spherical Radon 
transform of tp (see Appendix^ for the reconstruction ofTZ[ip] from its derivative). Note 
that the function tp itself can be reconstructed in a stable way from lZ[tp] using an inver- 
sion (filtered) retroprojection formula for the spherical Radon transform. From this, all 
inclusions are detected by the second statement in Lemma \2.b\ noticing that dAi is the set 
of discontinuous points of ip. 

6 A reconstruction algorithm of the true coefficient 

With all inclusions A\,A2,--- ,Ak in hand, we are able to find an initial guess for a* 
using the unique continuation property (Lemma I2.4|) and then employ a Landweber type 
iteration to reconstruct a*. As an initial guess, we reconstruct constant values inside 
each inclusion by minimizing the discrepancy between computed and measured boundary 
data. We prove a Lipschitz stability result for the reconstruction of the optical absorption 
coefficient in the class of piecewise constant distributions provided that the support of the 
inclusions is known. 



6.1 The data of boundary measurements and an initial guess 

Define 

S = j^c^l^ : a = a and , a*; € [a,a]| . 
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Let a\ and a 2 be in S. Their difference can be written as 

k 



a 2 -ai = y^/tjl y 



for some h = (hi, • • • , h^) G B = [a — a,a — a\ k . Note that B can be considered as a closed 
ball of M. k with respect to the oo— norm of M fc given by 

\h\ = max{|/ii|, • • • , 

The compactness of I? plays an important role in our analysis. Suppose that I ^ 0. 
Denote by $1 and $2 the optical energy density functions that correspond to a\ and a 2 . 
The function d> = $1 — $2 solves 



8=1 

+ = on <9f2. 



(6.1) 



Using as the test function in the variational form of (|6.ip . we see that 



(\V(j)\ 2 +a(j) 2 )dx + l / (d v <t>Yda< \h\X / |0|cfe, 
n Jan Jn 

where A is defined in Lemma 12.21 This implies 

\M\\v(aa) < C\h\ 

and, therefore, the continuity of the map h 1— > d u (p\gn. Since the map h G d^_k B 1 — y 
II^^IIl 2 (<9Q) is continuous and nonzero (due to Corollary 12. 5p . we can employ the com- 
pactness of d^kB in R fc to see that 

c ( a i) = m in > °- 

Identifying 5 with a compact subset of M fc , we can conclude that 



Properly scaling the inequality 



c = inf c(ai) > 0. 
ai€S 



for all /i G d^kB, we arrive at the following Lipschitz stability result using only one 
measurement. Note here that the support of the inclusions is known and only the value 
of the optical absorption coefficient inside each inclusion is to determine. 

Proposition 6.1 There exists c > such that for all a\,a 2 G Q, 

\\d v ®i - d„$> 2 \\L2(dn) > c\\a% - ci2||i,°°(n)> (6.2) 
where $1 and & 2 are the solutions of hl.l\) with a* replaced by a± and a 2 , respectively. 
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Remark 6.2 Inequality A6.2\) guarantees the uniqueness of the reconstruction for a* G S 
if ' d u $>*\oQ is considered as the data given. It, moreover, implies the stability in the sense 
that small noise does not cause large error. 

Proposition 16,11 suggests us to minimize the quadratic misfit functional: 

J(a) = ±\\d„$-d v $4l Hm , (6.3) 

where a varies in S and is the true optical energy density. This is possible since S is 
identical with a compact subset of M. k . By (|6.2p . the function aj = argmin J is close to a* 
provided that a* is a perturbation of a constant on each inclusion Aj. Therefore, aj can be 
considered as the background constant optical absorption distribution in the inclusions. 
For simplicity, we propose the following exhaustion method: for each fine partition P of 
the interval [a, a], try all values of on such that on equals each element of P, and finally 
choose the A:— tuple (ax, • • • , a^) that gives the smallest \\d u $ — d^^W^^y 

6.2 Internal data map and its differentiability 

Define the set that a* = (a*, • • • ,a* k ), identifying with true optical absorption coefficient 
a* of the form (jl.3p . belongs to 

k 

K := {a G JJ Wo' 4 (A,-) : a < a { < a and \\V ai\\ L 4, {A]) < 0,i = l, ■■■ ,k}, (6.4) 
i=i 

where 6 will be determined later in (|6.8p . It is obvious that K is closed and convex in H 
where # = U j= i H U A j) 

is a Hilbert space with the usual inner product 

(u,v)h = / Vnj • Vvjdx 
i=i Ja j 

for all u = (u\, • • • , Uk) and u = (ui, • • ■ , Vk) in H. 

Now, let the map F : K — > i7* be defined as follows. For all (oi, • • • , a^) G -fC, let 

a = ^ailA,, (6.5) 

i=0 

and 

F[a](u) = V / T[a] 2 V ai • for all v E H, (6.6) 

i=l 7 ' 

where T[a] was defined in Lemma |2.3[ We call F the internal data map. 
Theorem 6.3 T/ie map i 7 is Frechet differ entiable in K and 

DF[a](h,v) = ! (2T[a}DT[a}(h)Va i +T[a] 2 Vh i )Vv j dx (6.7) 
i=i jA i 
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for all a = (ai,-- ,«fc) £ K, h = (hi,-- - ,h k ) £ I\j=i ^^(Aj) n ^(Aj) and v = 
(v±, ■ ■ ■ , Vk) € H. Assume further 

0<9<^, (6.8) 

where ft' = Uj =1 Aj and Cqi is the norm of the embedding map of H 1 ^') into L 4 (f2'), 
multiplied with the constant in 12. 6\) . Then, DF[a] is well-defined on H and there exists 
a positive constant C such that for all h G H , 

\\DF[a](h)\\ H * > C\\h\\ H . (6.9) 



Remark 6.4 The term DT[a](h) in (6.1) is understood as DT[a] acting on the function 
that is equal to in Aq and to hj in Aj, j = 1, • • ■ , k. 



Proof of Theorem \6.3[ The Frechet differentiability of F and the expression (|6.7p of DF 
can be deduced from Lemma 12.31 and the standard rules in differentiation. We only prove 
IffTSl) . In fact, for all h€H, 

DF[a](h,h) = V f (T[a] 2 \Vhj\ 2 + 2T[a]DT[a](h)Va j Vh j ) dx 



>- E 

3=1 



{T[a] 2 \Vhj\ 2 )dx- [ \2T[a]DT[a}(h)VajVhj\da 



3=1 



> A 2 INIlr -El^ll^ r W( ft )ll i4 (^)ll Va jll i4 (^i)ll V/l ill i2 (^) 



It follows from the continuous embedding of H l (D) into L A (D) and (|2.6|) that 

DF[a](h, h) > A 2 (l - ||^ 1(D)5 



□ 



and therefore, inequality (|6.9p holds true. 

We now make use of Theorem 16.31 in order to prove a local Landweber condition which 
guarantees the convergence of the reconstruction algorithm. 
Let a and a' be in K. We can find t € [0, 1] such that 

\\F[a] - F[a']\\ H * = \\DF[ta + (1 - t)a'](a - a')\\ H * > C\\a - a'\\ H (6.10) 

by (|6.9p . Hence, if \\a — a'\\jj is small enough, then 

\\F[a] - F[a'} - DF[a](a - a')\\ H . < f]\\F[a] - F[a']\\ H * (6.11) 

for some r\ < ^. In other words, F satisfies the local Landweber condition (see |17j). 
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6.3 Landweber iteration 

Going back to equation (|5.ip . we have 

V • $ 2 Va = Aip (6.12) 

in the sense of distributions. However, the equation above can be understood in the classi- 
cal sense in each inclusion Aj. This observation plays an important role in reconstructing 
the true coefficient from the initial guess given in Subsection 16.11 
Considering Aip as an element of H * defined by 

- Aip(v) = ^2 f V ^ ' V V X > ( 6 - 13 ) 

for all v = (vx, • • • , Vk), we rewrite (|6.12p as 

F[a] = Aip. (6.14) 

Recalling that K is closed and convex in H, we can employ the classical Hilbert 
projection theorem to define the projection from H onto K as 

P : H 3 /i H> argmin{||/i - a\\ H : a G K}. (6.15) 

It is not hard to verify that 

\\P(h) - a\\ H < \\h - a\\ H (6.16) 

for all a € K. 

We next solve (|6.14p using the Landweber method to minimize 

I(a) = ±\\F[a]-Ai(>\\ 2 H ., 

where a varies in K with the initial guess ai = (a±, ■ ■ ■ , ctk), obtained in Subsection 16.11 
The corresponding guess for the coefficient is 

k 

ai = y~]ailAi- 

i=l 

There is a gap if we minimize I by the classical Landweber sequence given by 



a 



(n+i) = a (n) _ pDF[aW]*(F[a!ri] - Aip) 



because may not belong to K and Ffa^ 1 )] is not well-defined. Motivated by (|6.16p . 
which implies P (aW) is closer to o* than a^ n ' is, we modify this formula as 

fl (n+i) = p( a W) _ ^ J DF[P(a (n) )]*(F[P(a (n) )] - Atp). (6.17) 
We have the following convergence result. 
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Theorem 6.5 Suppose that the true optical distribution a* G K . Let be defined by 
{6.11 ) with a^ ) being the initial (piecewise constant) guess obtained as the minimizer of 
16. 3\) . Then the sequence converges in H to a* as n -f oo. 



Noting that F satisfies the local Landweber condition (see (|6.1ip ). we can repeat the 
proof of Proposition 2.2 in |17j to see that 

l|a (n+1) " + (1 " 2rj)\\F[P(a^) - A^]\\ 2 H , < \\P(a^) - o»|&. 

This and (|6.16|) imply 

||P(o(«+i)) _ a4 2 H - \\P(a^) - a4 2 H < (2 V - l)\\F[P(a^) - A^]f H , < 0. (6.18) 

It follows that 

oo 1 

£ \\F[P(aM)] - Wfn. < — — \\a4%, 
i=i ^ 

and hence 

F[P(a (n) )] -> Aip in J9#(fi) as n -> oo. (6.19) 

On the other hand, we can see from (|6. 18|) that the sequence (P(a( n ))) n >i is bounded 
in H. Assume that P(a^) converges weakly to a' for some a' G H. Since K is closed and 
convex, it is weakly closed and therefore a' G K. Passing to a subsequence if necessary, 
this sequence converges to a' a.e. and also converges strongly to a' in n^=i L 2 (Aj). So, 
T[P(a n )} converges to T[a'] in H 1 ^) and hence in L (Q). For all v G H, we have 

f (T[P{a^)} 2 VP(a {n) ) -T[a'] 2 Va')Vvdx 

3=1 JA = 
k 

= J2 I (T[P(a^)} 2 -T[a'} 2 )VP(a^)Vvdx + f T[a'] 2 {VP{a^) -Va')Vvdx 
j=l IJAj J Aj 

which goes to by the dominated convergence theorem and the weak convergence of 
P(a ( - n ^) to a' in H. We have obtained F[a'] = Aif) = F[a*]. Using (|6.9j) gives a' = a*. 

In summary, if the true coefficient a* is a perturbation of a constant on each inclusion 
then the coefficient ai obtained in Section [6.11 is quite closed to a*. Moreover, the misfit 
between the initial guess aj and the true distribution a* can be properly corrected by the 
sequence in (|6.17p . 



7 Concluding remarks 

In this paper we have introduced a Landweber scheme for reconstructing piecewise smooth 
optical absorption distributions from opto-acoustic measurements and proved its conver- 
gence. Because of the jumps in the absorption coefficient, we have used weak formulations 
for the Helmholtz decomposition for & 2 Va* and the relation between the spherical Radon 
transform of its gradient part V an d the cross-correlation measurements M v (y,r). Note 
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that we can enrich the set of data as follows. For / € L 2 (<9S7) such that / > a.e. on <9S7, 
compute instead of (|3.4|) the quantity 

Mf*{y, r) = \ [ (Jdv&n - gd v &) da, y € S^r > 0. 
V Jdn v ' r 

Similarly to f)3.5|) . integration by parts yields 




where is the solution of (jl.ip with 5 replaced by /. 

The enriched data (|7.ip may be used in order to generalize our approach to the case 
of measurements of the outgoing light intensities on only part of d£l by choosing / sup- 
ported only on the accessible part of the boundary. Another interesting and challenging 
problem is to prove statistical stability of the proposed reconstruction with respect to 
a measurement noise by combining Fourier techniques together with statistical tools [8]. 
Numerical implementation of the Landweber-type iteration is under consideration and will 
be the subject of a forthcoming publication. The behavior of the proposed method with 
respect to the optical absorption contrast will be investigated. It is expected that higher 
the contrast, more efficient the method is. 

A Proof of Lemma 12.41 

The boundedness of (j> together with the assumption that = on dtl imply by standard 
regularity results that <p € C 1 (9f2 U f2 \ ft ). Arguing similarly to Proposition 2.5 in [5], we 
see that G C 2 (fl \Tt'). Define 

U = {x en\Tl' : u{x) + 0}. 

The continuity of 4> shows that U is open. Assume, on contrary, that U is nonempty. 

Noting that U can be decomposed as the union of its connected open subsets. Denote 
by O the connected component of U, which is closest to d£l. Without loss of generality, 
assume that > in O. Let 

5 = dist(0,dO). 

The distance above is understood as the length of the shortest curve, contained in f2 \ fi 
and connecting O and dQ.. 

In the case that 5 = 0, dO and dQ have a common point xq. Applying the Hopf lemma 
for the equation 

-Acf) + c(f) = in O, 
<t> > on dO, 

gives d u (p(xo) < 0, which is impossible. 

When 5 > 0, it is easy to see that 4> = in a neighbourhood of dVL. Assume that such 
a neighbourhood and O have a common boundary point xq. Noting that V0(xo) = 0, we 
can apply the Hopf lemma again to get the contradiction. □ 
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B Proof of Proposition 



32 



We write u and v when referring to Uy jT and Vy )T respectively for simplicity. Using (|3,7|) 
and the same arguments when estimating \\a u — a||x,i(n) m the previous section yields 

IK - a||L2(n) < 0(rj). 

This, together with standard i? 2 -regularity results (see, for instance, |16[ Theorems 8.8 
and 8.12]) and the embedding of H 2 (Vl) into L°°(fl), gives 

||*«-*IU«(n) < 0(rj). 

Hence, it follows from (13.61) that 



a)&& u dx - / (a u — a)<& 2 dx 


< 


J 


Jn 




Jn 




< 


11$ 



< ll^llL°°(n)ll a M - allien) II 3>« - *||ioo(n) 



< erf 



(B.l) 



The constant c depends only on a = max a and a = mina, both of which are assumed to 
be known. The independence of c on ||$||i,°°(f}) can be deduced from Lemma 12.21 Now, 
note that the second integral in the left hand side of (|B. 1|) can be rewritten as 



/ (a u - a)<& 2 dx = ( a u — a)<& 2 dx 

Jn JAiUP(Ai) 



E 

i=l L 



/ (a u - a)§ 2 dx + (a u - a)$ 2 dx 

JA i C\P{A i ) JAiAPiAi) 



AiAP(Ai) 



and that the integral in (|4.ip is equal to 

f(a- a )V • (§ 2 v)dx = V / (a- a )V • ($ 2 u)da 
Jn ~^ J a. 



y~ y j / (aj - a Q )$\ ■ Vida - / $ 2 Va • vdx 
i=1 UdAi J Ai 

n r- 

^ (ai- a )$ 2 v ■ Vido 
~T YJdAi 

$ 2 Va • vdx 



AiC\P{Ai) 



L 



$ 2 Va • vdx 



Ai\P(Ai) 



Therefore, we have 

/ (a u - a)$ 2 dx + (a- a )V ■ (§ 2 v)dx 
Jn Jn 

j (a u - a + Va • u)$ 2 dx + / $ 2 Va* • 

iA^PCAi) JAi\P(Ai) 



i=l 



+ 



A,AP(4,) 



a u — a)$ 2 <ix — / (ai — ao)$ 2 f ■ Vida 
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Denote by a«, ft and ji the last three quantities in the inequality above. We need to prove 
that they all are bounded by 0{rf) to complete the proof. 

(i) Since for all i = 1, • • • ,k, a,i G C 2 (Ai) and H-D^H^co^.) are bounded by some known 
constants, we can find a constant c such that 



(a u — a — Va • u)<& 2 dx 



'AinP(Ai) 

On the other hand, using the classical substitution method in integration gives 



< 



Va • (u + v)$ 2 dx 

! ($ 2 Va)(£)- I'" (l Ai (u + v))((l + ^)dpd^ 
Js(y,r) J-v r 



+ 



S(y,r) J-T] 



($ 2 Va) + (cD 2 Va)(0 . + v)) (1 + ^ dpd£ 



<0 + 2r7 2 ||d r ($ 2 Va)|| L o O(A) |£ ?7 | 
< cr/ 3 . 

The quantity 04 is bounded from above by 0(rj 3 ) because 



on < 



/ (a u - a - Va • u)<£> 2 dx + / Va • (u + v)$ 2 dx 



(ii) The fact that ft < 0(rj 3 ) can be deduced from the boundedness of the integrand 
and (1X711 . 



(in) The main point of the proof is the estimate of 7$. Denote = f when 
it is defined and 

X-.dAx [0, 77] — ► AjAP(Aj) 
(z, i) 1 — >• z + tv(z). 

X is well defined when ■y(z) is non-zero and not parallel to dA{. For any z £ <9.4j 
satisfying this condition, we denote T(z) the tangent plane to dAi and -B(z) a basis 
adapted to the sum R d = T(z)®Rv(z), Let I^-i denote the (d — 1) x (d — 1) identity 
matrix. We get 

"J d _l + idu(z) 0" 
* 1 



dx(z,t) 



and as u (z) = (z — y)/\z — y\. The operator dv{z) does not depend on r\ and tdv(z) 
0{rf) with a constant depending on tq. Then, 



det(dx(z, t)) = l + tV ■ (v)(z) + 0(7? 2 ) = 1 + Q(rf). 
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As B(z) is not orthonormal, the differential volume written with the variables (z, t) 
depends on the angle between v{z) and v(z) called 9{z). This volume at the point 
z + tv(z) is (1 + 0{rj)) cos{9{z))dtdz. Knowing this, we denote 

(dA^ = {z£ dA h ±9{z) > 0} 



and write 

r r rH z )\ 

/ {a u -a)<S> 2 dx = / / (a u -a )<S> 2 (z + tv(z))(l + O(ri)) 

JP{Ai)\A t J{dA,)+ JO 

x cos(6(z))dtdz 

and as a\ and 3> are C 1 (Aj), we can write that for any z € (dAi) + , and t € [0, \v (z)\], 
\(a u - a )$ 2 (z + tv(z)) - (a* - a )<S> 2 {z)\ < 0{rj). 

Then, 



v(z)\ 



(a u -a )<5> 2 (z + tv(z))(l + 0(r]))dt - (m - a )<S> 2 (z)\v(z)\ 



< 0(r? 2 ). 



Now, noticing that cos(9)\v(z)\ = v(z) ■ u(z) and that a((dAi) + n S^), the surface measure 
of (cMj) + n E„, is of order 0{rj), we have 



We also get 



/ (a u - a)$ 2 - / (at - a )^ 2 v ■ v 

JP(Ai)\Ai J{dA l )+ 



/ (a u - a)$ 2 - / (oj - a )$ 2 w • i/ 



< 0(r? 3 )- 



< C>(r/ 3 



by the same arguments. 



□ 



C Construction of 1Z[ip] from formula (15.31) 



In order to construct from formula (|5.3p . we need to invert the operator : 

L 2 (C) — > G~ 1 (C) and prove the stability of the inversion. For any / G L 2 (C), by 
Fubini's theorem, the function F(y, r) = f Q r f(y, p)dp is well-defined and in G(C) but not 
in Gq{C). Since this operator is acting on distributions which are zero on S^xJO, tq[, we 
introduce 



P : L 2 (C) 



G (C) 
(y,r) H' 



R 

v(y,p) - — X]o,r {(p)v(y>P R / r ^ ) d P 



and its dual 



p* : G -1 (C) — ► L Z (C). 



(C.l) 



The following result holds. 



22 



Proposition C.l For all f £ L 2 (C) such that f = on S^x]0, tq[, we have the inversion 



formula 



Proof. For any ip £ L 2 (C), we have JfMv 9 ] = — on x [ro, R[ and therefore, 



2i 

dr 



G-i(C),Gj(C) 









/ p* [«] p 













which yields the claimed result. 

Proposition C.2 For all u £ M := G G _1 (C) : supp(u) C S^x [r ,R[}, 

\\p*u\\mo < \WWg-hc) 

Proof. We first note that p* [u] = on 

SnX]0,r o [. Then, for any tp £ L 2 (C), wc set 
P*MX[r ,R[V < ll«llG- 1 (C0l|P[X[r O , J B[¥']ll(3j(C) 

a . 

^ H u llG-i(G)l!^[X[ro,R[^llL2(G) 

< hllG-i^llXtro^lU^C) 

< IMIg-hcjIMIl^c), 

and the proof is complete. 

Finally, we deduce the following result. 

Corollary C.3 From formula \5.3\) . we have 

1 



□ 



□ 



p*(r d - 2 M). 



Moreover, for n small, if TZ[ifj v ] 



ro\\w\\l 

1 „*f^d-2 



ro \\w\ 



»-^]||^(C) 



-p*(r d - 2 M v ), then 

R d-2 



< 



M - M v \\ G -i (c) , 



ro\\w\\i 

which insures the stability of the construction oflZ[ip] from the measurements M„. 
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